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NOMENCLATURE 


r,e,z  cylindrical  coordinates  used  to  describe  the 

undeformed  configuration  of  the  plate 

h  thickness  of  the  plate 

a  outer  radius  of  the  circular  plate 

b  outer  radius  of  the  mass 

p  mass  density  of  the  plate 

t  time  variable 

u,w  radial  and  transverse  displacements  of  the 

middle  plane,  respectively 

€    ,e„  total  strains  in  the  radial  and 

r  e 

circumferential  directions,  respectively. 

o    o 

€  ,€^  radial  and  circumferential  strains  acting  on 

r  e 

the  middle  plane,  respectively 

c    ,c^,c  radial,  circumferential  and  normal  stresses, 

r   e   z 

respectively 
N  ,N  membrane  forces  per  unit  length 

T  ,T  kinetic    energy   of    the    plate   and   of 

p  m 

concentric  rigid  mass,  respectively. 

U  , U,  strain   energy  due   to   stretching   of   the 

middle  plane   and  due   to  bending  of   the 
plate,  respectively 

W  work  done  on  the  plate  by  external  forces 

V/,^)  stress  functions 

M  ,M  bending  moments  per  unit  length 

r   0 

q(r,t)         time-variant  loading  intensity 


Iv 


E  ,E  ,E        material    constants    in   the   directions 
r  e  re 

indicated  by  subscripts, 
a   ,a   ,a      elastic  constants  (compliance  coef f iciences) 

R  radius  ratio,  =  b/a 

y  mass  ratio 

a  elastic  constant  ratio,  =  E^./E 

r6   r 

B  elasitc  constant  ratio,  =  E./E 

C  elastic  constant  ratio,  =  E  /E 

r   " 

V  Poisson's  ratio,  =  a/fi 

3 

D  flexural  rigidity  of  the  plate,  =  Eh  /12 

I  action  integral  for  the  vibrating  system 

SI  first  variation  of  I 

n,C  admissable    variations   of   w   and   u, 

respectively 

€  arbitrary  infinitesimal  constant 

^,T  dimensionless  space  and  time  variables, 

respectivley 

5W  virtual  work  of  transverse  forces 

X  dimensionless  transverse  displacement 

Q(?).Q  (?)  dimensionless  loading  densities 

9{?)»f(?)  shape  functions 

A,Af  amplitude  parameters 

X  nondimensional  eigenvalue 

u  nondimensional  angular  frequency 

cj,  linear  angular  frequency 


V 


y,Z,H  (6x1)  vector  functions 

[M],[N]         coefficient  matrices 

"ri  missing  initial  values 


VI 


Chapter  I   INTRODUCTION 

Anisotropic  materials  play  an  important  role  in  modern 
technology,  because  lots  of  new  anisotropic  materials,  such 
as  reinforced  plastics  and  composite  materials,  are  being 
used  in  missiles,  aircraft,  space  vehicles,  pressure 
vessels  and  parts  of  structures  to  meet  the  special 
requirements . 

In  the  past,  materials,  regardless  of  their 
composition,  were  usually  considered  to  be  homogeneous  and 
isotropic  because  such  assumptions  make  calculation  simple. 
However,  these  simplified  assumptions  often  lead  to 
inadequate  or  incorrect  results.  Today,  to  meet  the 
sophisticated  technology  requiremehts ,  we  need  to  consider 
the  anisotropic  property  of  materials,  that  is,  the 
differences  in  elastic  properties  of  materials  in  various 
directions. 

A  circular  anisotropic  plate  whose  material  symmetry 
of  aeolotropy  at  a  point  is  in  radial  and  tangential 
directions  is  called  an  orthotropic  plate. 

The  governing  differential  equations  describing  the 
motion  of  thin  circular  plates,  in  general,  are  nonlinear 
and  coupled.  Due  to  the  complex  nature  of  the  resulting 
governing   equations,   an  analytical   solution    is    very 


difficult  to  obatain.  Thus,  approximate  methods  [1]  of 
analysis  must  be  used. 

Large  amplitude  vibrations  of  a  clamped  isotropic 
circular  plate  with  a  concentric  rigid  mass  at  the  center 
have  been  studied  by  Handelman  and  Cohen  [2],  Chiang  and 
Chen  [3]  and  Becker  [4]. 

On  the  other  hand,  large  amplitude  vibrations  of  a 
clamped  orthotropic  circular  plate  with  isotropic  core  have 
been  studied  by  Huang  [5].  Woo  [6]  payed  attention  to  the 
same  problem  except  he  limited  the  vibration  to  a  small 
amplitude.  However,  the  addition  of  a  concentric  rigid 
mass  attached  rigidly  to  the  orthotropic  circular  plate  has 
received  limited  attention. 

In  this  thesis  Hamilton's  principle  is  applied  to 
derive  the  basic  differential  equations  for  the  problem  of 
the  large  amplitude  of  an  orthotropic  circular  plate 
carrying  a  concentric  rigid  mass  at  the  center.  Assuming 
the  system  performs  the  harmonic  vibrations,  the  time 
variable  is  eliminated  by  employing  the  Ritz-Kantorovich 
averaging  method.  Therefore,  the  basic  governing 
differential  equations  are  reduced  to  a  pair  of  ordinary 
differential  equations,  which  form  a  nonlinear  boundary- 
value  problem  with  the  associated  boundary  conditions.  For 
the  purpose  of  numerical  computation,  this  boundary-value 
problem   is  converted  to  the  related  initial-value  problem. 


The  relation  among  natural  frequency,  stiffness  parameters 
and  mass  ratio  for  free  vibration  as  well  as  for  forced 
vibration  are  investigated,  and  presented  by  figures. 
These  results  provide  useful  inf  Oiir.ati  on  for  the  design  of 
the  nonlinear  vibration  of  an  orthotropic  circular  plate 
with  a  concentric  rigid  mass.  The  free  and  forced 
vibration  with  variable  parameters  are  also  investigated; 
their  effects  and  the  corresponding  results  are 
illustrated.  Finally,  the  radial  and  circumferential 
stresses  in  free  vibration  are  also  determined.  These 
results  of  stress  distribution  also  provide  useful 
information  for  stress  analysis  in  industrial  utilization. 


Chapter  II    DERIVATION  OF  GOVERNING  EQUATIONS  OF  MOTION 

Consider   the   thin   circular  plate  with  a  concentric 
rigid  mass  M  as  shown  in  Fig.   1.   The  plate,  except  for  a 

concentric   rigid  mass  M   of  radius  b  occupying  the  central 

portion,  is  composed  of  an  elastic,  homogeneous  and 
cylindrically  orthotropic  medium  bounded  by  the  planes 
2=+h/2.  The  plate  is  clamped  at  the  outer  edge  of  radius 
a.  In  discussing  circular  boundaries  problems,  it  is  more 
convenient  to  use  cylindrical  coordinates  (r,  e,  z)  and  set 
the  origin  at  the  center  of  the  middle  plane  to  describe 
the  motion. 

Before  deriving  governing  differential  equations,  we 
have  made  the  following  assumptions 

1.  Any  straight  line  normal  to  the  middle  plane  before 
deformation  remains  a  straight  line  normal  to  the 
neutral  plane  during  deformation. 

2.  Normal  stress,  o      is   small   in  comparison  with  the 

other  stresses,  and  may  be  neglected. 

3.  Circular  symmetry  with  respect  to  the  z-axis  is 
retained  during  vibration. 

4.  The  effects  of  stretching  the  middle  plane  of  the 
plate  are  included. 


orthotropic  plate 
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Fig.  1.  Circular  Plate  and  Cylindrical  Coordinate  System 


5.  The  vibration  of  the   plate   is   not   limited   to   the 
infinitesimal . 

The  radial  and  circumferential  strain-displacement 
relations  derived  from  large  deflection  theory  [8,10]  are 
as  follows: 

1   2 
e=u  +  -w-zw  (la) 

r   r   2   r     rr 

e^=  H  _  z  w  (lb) 

e  r   r  r 

where  u(r,t)  and  w(r,t)  denote  the  radial  and  transverse 
components  of  the  mid-plane  displacement,  respectively, 
while   the  subscripted  variables  w  and  w   represent  the 

first  and  second  partial  derivatives  of  w  with  respect  to 
r . 

The  strain-stress  relations  in  the  polar   coordinate 
system  are  : 


cf«=  E^€^  +  E  ^€  (2a) 

e   e  e    re  r 

a  =  E  ^€^  +  E  €  (2b) 

r   re  e    r  r 

where  <t  and  cr  are  normal  stresses  in  the   tangential   and 

radial   directions   respectively,   E„ ,   E  ^   and  E   are  the 

^      -^    e    re      r 

material  constants  (Young's  modulus),  and   e    and   €   are 


normal   strains   in   the   tangential   and  radial  directions 
respectively. 

The   stress-strain  relations  in  equations  (2)  can  also 
be  represented  in  the  form  [9] 


S=  ^ll'^e  -^  ^12^r 


'r=   ^12*^6  ^  ^22^r 


(3a) 
(3b) 


where  a. .  are  elastic  constants. 

Equations  (3a)  and  (3b)  can  be  written  in  matrix  form 


11    12 
^12   ^22  J 


r  € 


e 


Let 


A  = 


11     12 
L  ^12    ^22 


X  = 


e 


B  = 


the  matrix  can  be  expressed  symbolically 

AX  =  B. 
then 


X  =  A~^  B 


which  implies 


r  J 


22 


^11^22  ~  ^12 


-  a 


12 


-  a 


12 


^11^22  "  ^12 


^11^22  "  ^12 


11 


^11^22  "  ^12 


thus  a      and  a      can  be  represented  in  terms  of  a.  .  ,  €   and 


6   by 


^e  = 


22 


e_  - 


12 


2    ^e  2    r 

^11^22    ^12         ^11^22  ~  ^12 


(4a) 


-  a 


a   = 

r 


12 


6_  + 


11 


2         "e    '  2    r 

^11^22  "  ^12         ^11^22  "  ^12 


(4b) 


and 


^e  = 


22 


^11^22  "  ^12 


(5a) 


E   = 
r 


11 


^11^22    ^12 


(5b) 


-  a 


12 


re  2 

^11^22"  ^12 


(5c) 


Here  we  define 


a   = 


^re  ^12 


11 


$  = 


'e  "22 


11 


c    = 


E  a,  ,  , 

rill 


^e  ^22  ^ 


1/  =  - 


12 


22 


/8 


Integrating  the  in-plane  stresses  across  the  thickness 

of  the  plate,  radial  and  circumferential  bending  moment  per 
unit  length,  M  and  M  ,  respectively,  are  found  to  be 

r    J-h/2   r 

=  f"^J?^J  (E  €   +  E  ^6^)zd2 
J-h/2    r  r     re  e 

=    I'^m    [E    (u   +  -^  w^    -    zw      )+E   ^(-   -   -  w    )]zdz 

J-h/2    "■    r^    r  2      r  rr        re'r        r      r    ■" 

^3 

= v^r-    (E   w  +   -   E    ^w    ) 

12  r    rr        r      re    r 

=   -   D    (w        +  -  aw    )  {6a) 

rr        r        r 


f+h/2  , 

^  =  J-h/2  ^e  ^^^ 


=    I'^^il    (E    ^6      +    E^6^)zdz 
J-h/2    ^    re  r         6  e 


=    f'^^^'/o    [E    „(u   +   ^  w^    -    zw      )+E^(-   -   -  w    )]zdz 
J-h/2    "-re^r      2      r  rr        er        r      r    ■' 

^^'    r  ^re         ^  ^e,^r,i 

= 12-    L-E—  "rr-^   r^~U 

r  r 

=   -   D    (a   w        +  i   /S  w    )  (6b) 

rr        r  r 

3 
in  which,  D=E  h  /12  is  the  flexural  rigidity  of  the  plate. 

In  addition,  a     and  /9   has   the   following  restrictions 
[Appendix  A]  : 

0  <  a  <  1 


and  a  <  /S 

In  the  same  way,  the  expressions  for   the   radial   and 
circumferential   forces   per   unit   length,   N^  and  N^ ,  are 

obtained  by  integrating  the  respective  stresses  across   the 
thickness  of  the  plate, 

T        J-h/2   r 

=  I-h/2  (^r^  ^  ^re^e^  ^- 
=  E^(u^h  +  I  w2h)+E^^(fh) 

=  iiD       1   2     u  (,^) 

,  2    r   2   r     r' 
n 


=  E^e^V^  i  w2h)+E^{Hh) 

=  ^   [a(u  +  \   w2)+^  H],  (7b) 

,2     r   2   r    r 
n 


The  Energy  Method  : 

To  formulate  the  boundary- value  problem  ,  we  shall  use 
Hamilton's  principle.  The  power  of  this  approach  is  that 
it  furnishes  automatically  the  correct  number  of  boundary 
conditions  and  the  differential  equation.   In   this   system 


10 


during   the  motion,  three  kinds  of  energy  and  the  work  done 
by  the  external  forces  are  considered  as  follows: 
U, =  strain  energy  due  to  bending 

U  =  strain  energy  due  to  stretching  of  the  middle  plane 

T  =  kinetic  energy 

W  =  work  done  by  the  external  forces . 

The   strain   energy  due  to  bending  can  be  expressed  as 
[10] 


U.=   -fj'^f?    (J  M  w      +  -|-  M^w    )rdrde 
b        J0Jb^2      rrr        2r      er 


Substituting  equations  (6a)  and  (6b)   into   this   equation, 
one  obtains 

U,=  -x    L.  L.[(w  +  —  w  )w  +(aw  +— w  -w  ]rdrde 
b   2  Jo  Jb"-   rr   r   r'  rr    rr   r  r  r  r' 

=  rrD  f^  (w  ^+  2-w  w  +  ^„w^ )  rdr  (8) 

J  b   rr   r  rr  r    2  r 

r 

Strain  energy  due  to  stretching  of   the  middle  plane 

can  be  expressed  as 


U   =    \V^\^    (^N   €%  ^N^€*)rdrd0 
s      J0Jb^2rr      206 


substituting   equations   (7a)   and  (7b)  into  this  equation, 
one  obtains 


fa  fl2D, ,   ,  1  2 ,  u. ,   .12,, 

U  =  T  ,   -^ [  ( U  +  -W  +a- )  ( U  +7rW  )  ] 

s    Jb  I  .2'-*  r   2  r   r' ^  r  2  r'-" 
h 

lcL^^^^\^l)^^{l)Ul)]   rdr 


h 
12D 


11 


=  ^[^[u^u  w^sju   uAw^^  n^iw^rdr  (9) 

,2   Jb      rrrrrrrz        4r 
n  r 

o  o 

where  e   and  e„  represent  the   radial   and   circumferential 
r      e   ^ 

strain  at  the  mid-plane,  respectively. 

Kinetic  energy  due  to  translation  can  be  divided   into 

two   portions,  one  is  the  contribution  of  the  plate  portion 

T  ,  the  other  one  is   due   to   the   concentric   rigid   mass 
P 

portion  T  .    The   radial   component   of  velocity,  in  this 

problem,  is  rather  small  in  comparison  with  the  transverse 
component  of  velocity  and  can  be  neglected.  Thus,  only  the 
transverse  component  of  velocity  is  considered  in  kinetic 
energy. 

T  =  Iq'^JJ  -(phw^)rdrde  =  nph   J^  w^  rdr  (lOa) 

where   p   is  the  mass  per  unit  volume  of  the  plate,  and  the 

kinetic  energy  T   is 
^'   m 


1     2 

T  =  4  M  w^ 

m   2   c  t 


,  (10b) 

r=b  ' 


where  M   is  the  concentric  rigid  mass, 
c  ^ 

The  work  done  by  the  external  forces  can  be  expressed 

as 

W  =  Jo"'JJ  q(r,t)w{r,t)rdrde  (11) 

where  q(r,t)  is  the  time-dependent  external  loading   inten- 
sity, assumed  to  be  symmetric  with  respect  to  the  z-axis . 


12 


Consider  the  motion  of  the  plate   system   between   two 
fixed  instants  t   and  t  ,  and  applying  Hamilton's  principle 

,  the  first  variation  of  the  action  integral  61  must  be 
zero . 

I  is  the  action  integral  defined  by 

1=  111    (T  +T  -U^-U,  +W)dt  (12) 

J  t-   p  m  s   b 

substituting  equations  (8),  (9),  (10a),  (10b)  and  (11) 
into  equation  (12)  and  taking  the  first  variation  of  I,  one 
obtains  61. 


61=   -27rpheJ^lJ^  )7w^^rdrdt-€M^J*l  rjw^^dt 

-2rrD€j^l  n^.  (  rw^j.+orw^  )  |  Jdt  +  2;rD€j^l  H  (  rw^j.^+w^^ )  |  ^dt 

-27rD€jJl  '?{fw^)|bdt  -  2;rD€j^lj^  n  ( rw^j.j.j.+2w^^^)drdt 

+  2rrD€f^lf^  r\{^w      -   ^w  )drdt 
J  t^J  b  '  r  rr    2  r 


"0 

12n-De    ft 


h^      -0 


J^l   C(2ru^+rwJ+2au) l^dt 


12^06  ft 


h^    -0 


ft  3  la 

Ll    n(2ru  w  +rw   +2auw    )Ldt 
J  t^    '  ^        r   r        r  rib 


+    ^2n-Dert    ra   2^{ru     +u  +  ^w^+rw  w     +au    )drdt 
.2   Jt^Jb     ^^      rr      r      2   r        r   rr        r 
h  0 


127rDe  ft,  fa 


j^lj^  C(2au^+aw^+|/8u)drdt 


h^   '  -O'b 


13 


127rD€  ft,  fa 


[^lj^)7l2w    [u    (l+a)+ru   ^]+w^^(  2au+2ru^+3rw^  ) +w^|drdt 


•Ob 


+  27reJ^lJ^  rjq(r,t)rdrdt 
0 


=  0 


The  necessary  condition  to  let  <5I=0  is  that  the  in- 
tegrands of  the  double  integrals  and  the  single  integrals 
must  vanish  separately.  Thus,  the  double  integrals  provide 
the  governing  differential  equations  of  motion, 


D(w    +  -w   -  -^w   +  -^w  )+phw.  . -q(r,  t) 
V  rrrr   r  rrr    2  rr    3  r   "^   tt  ^ 

r       r 

=  — ;r(— u  w  +u  w      +u     w  +-U  w  +-UW      +  77— w   +  — w  w      )  (13a) 

2^r   rr      rrr      rrrrrrr      rr      2rr      2r   rr 


and 


112        a      2      B 

U   +— U  +-T—Vi    +W  W   -  -r—W    -    —rU    =     0 

rr  r  r  2r  r   r  rr   2r  r    2 

r 


(13b) 


while  the  single  integrals  yield  the  boundary  conditions 


27rDr(w      +— w    )n 
rr   r   r    'r 


=   0    , 


247rD       ,       ^1    2    a    ,_ 

r  (  U    +-;rW    +— U    C 

2         ^    r    2    r    r    '^ 


=    0 


(14a) 
(14b) 


1  fi 

2/rDr{w        +   — w      -   — ^w    )n 
rrr      r    rr        2    r 
r 


a         24n-D         ,         a      1    2.       a 

.  -   ;;— rw    (u    +— U+-T-W    )n    , 

b  ^2        rrr      2r'lb 


-    M^w^t^b   =    ° 


(14c) 
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The  corresponding  geometric  boundary  conditions  for  a 
clamped  orthotropic  circular  plate  carrying  a  concentric 
rigid  mass  are 


w    =0        w.=  0       w     =0 

r=a  r I r=b  rlr=a 


u    =  0         u   ,  =  0  (15a) 

' r=a  I r=b 


which  means  transverse  displacement  at  r=a  and  slope  both 
at  r=a  and  r=b  are  equal  to  zero,  and  radial  displacements 
at  mid-plane  are  also  zero  at  r=a  and  r=b. 

The  geometric  and  natural  boundary  conditions  are 
supplemental  to  each  other  and  add  up  to  the  right  number 
of  boundary  conditions.  Their  satisfaction  ensures  that 
the  solution  of  the  differential  equation  is  unique. 


Since  W    ,.=  0  and  W 
r I r=b         r 


=  0,  equation  (14c)   can  be 
r=a      ^ 


reduced   to  be 


1  ^ 

(w        +  —  w      ) 
rrr      r      rr 


■^        w^^L  _w  (15b) 


b  27rbD        ttlr=b 

Now,  introduce  the  stress  function  \/^(r,t) 

12D      r  2 

>fi{T.t)    =  2^^"r''"  2'^r.'''"^^ 

h 

which  has  the  following  relations  with  N^  and  N. 

^  =  N  I*  =  N^  (16a) 

r    r  9r    © 

where  N  and  N  must  satisfy  the  equilibrium  equation  [11] 
r      0 
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aN       N  -N. 
-X-^-   +  —^ ^  =  0  (16b) 

Using  these  relations,  equation  (13a)  and  (13b)  can  be 
transformed  as  follows  : 


D  w    +  -w   -  -^w  +  -f  w   +Phw^^-q(r,t   =  -  w  \y)  17a) 

V  rrrr   r  rrr    2  rr    3   r'  '^   tt  ^^    '    r   r^'r   ^    ' 
r      r 

^  +   k^    -   L^   ^    = J:_  E  h(/9  -  a^)w^  (17b) 

r 


At   the  mid-plane   of   the  plate,  the  circumferential 

o  o 

strain,  €„  can  be  expressed  as  €.=  u/r.  Then  from  equations 

o  © 

(2),  (7)  and  (16a)  we  can  obtain 

rE 

u  =  ^-y—  (v//  -  fv^)  =  0 

E  E«-E^ 

re  re 
The  conditions  for  u  to  be  vanished  at  r=a  and  r=b  are 

(,/, ^t)   .  =  0 

For  the  purpose  of  simplification  in  mathematics,   the 
following  dimensionless  quantities  are  introduced 

2 

X   =   w/a    ,  ?    =   r/a    ,         R  =   b/a    ,         y   =   M^/7rb   ph    , 


T   =   t(/SD/pha^)^,  ^   =   a^g^/ha    , 
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and 


Q  =  12a22(|)'q 


Then,  the  governing  differential  equations  (17a)  and  (17b) 
and  boundary  conditions  (15a)  and  (15b)  are  converted  into 
the  following  non-dimensional  forms  : 

2 


^????   ?^???~   ^2^??' 


73^?^  ^x^^-(i-  |-)Q 


f(i)^l- 


-j-)ix^4>)^    . 


*??^  f*r  -7^  *  -ffx^  =  0  , 


(18a) 
(18b) 


and     (  X  )t=i  =  ° 


<  *?  -  ?<^  )?  =  1  =  ° 


(19) 


(  X^^^+  |X^^+  2^X^t)?=R  =  ° 
Equations  (18)  and  (19)  form  a  set  of  nonlinear  bound- 
ary value  problems  describing  the  large  oscillation   of   an 
orthotropic   circular  plate  with  a  concentric  rigid  mass  at 
the  center  of  the  plate. 
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Chapter  III   APPROXIMATE  ANALYSIS 

The  differential  equations  of  motion  together  with  the 
associated  boundary  conditions  constitute  a  boundary-value 
problem.  Any  solution  to  the  problem  defined  by  the 
nonlinear  differential  equations  must  satisfy  the  boundary 
conditions  .  It  is  usually  very  difficult  to  obtain 
analytic  (closed-form)  solutions  to  the  boundary  value 
problems,  thus,  the  approximate  method  is  usually  applied. 
The  Ritz-Kantorovich  averaging  method  [7],  which  has  been 
successfully  applied  to  numerous  elasticity  problems,  is 
adopted  to  eliminate  the  time  variable  and  reduce  the 
governing  differential  equations  to  a  system  of  nonlinear 
ordinary  differential  equations,  which  are  the  only 
function  of  the  space  coordinate. 

The  Ritz-Kantorovich  Method 

Assuming  the  dynamic  system  performs  the  harmonic 
vibrations,  then  the  functions  Q(?  ,t)  ,  x  (? -t)  and  <|>(?,t)  of 
equations  (18a)  and  (18b)  can  be  expressed  as 

Q(?,T)=  Q*{?)(5)sinuT  (20a) 

X(?,T)=  A  g(?)sin(JT  (20b) 

<^(?,T)=  A^f(^)sin^uT  (20c) 
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where  A  denotes  a  nondimensional  amplitude  parameter;  u      is 
the   nondimensional   angular   frequency;   g(^)  and  f{^)  are 

shape  functions   to   be   determined   corresponding   to   the 

* 
functions   x   and  <p ,      respectively,   and  Q  (^)   is   the 

dimensionless  loading  function.  Since  equations  (20a), (20b) 

and  (20c)  can  only  satisfy  equation  (18b)  ,  but  not  satisfy 

equation  (18a)  exactly  for  all  t,  thus  the  residual  may  be 

found  and  minimized  by  the  use  of  Ritz-Kantorovich  method. 

For  any  instant  of  the  dimensionless  time  t,   the  virtual 

work   of   all   the  transverse  forces  as  they  move  through  a 

virtual  displacement,  <5x  =  AsincjT  ( <5g)  ,  is 

dW  =  J^  k^6x    ?  d?  (21) 

where 

*i  =  "nf?-^  ^nr  -fjXf?^  -^x?*  *'(rr-<i-  ^i" 

It  is  reasonable  to  require  that  the  integral  of  this 
virtual  work  over  a  complete  period  of  vibration  be  zero 
[12] ,  that  is, 

Jo"'''"  (<SW)dT  =  0  .  (22) 
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Substitution  of  equations  (20a),  (20b),  and  (20c)  into 

* 
(21),  (5W  can  be  represented  in  terms  ofA,  g(^),   f(?),   Q 

(?)  and  sinuT, 


dW=  sin2uT{J^[g^^^^.|g^^^-  ^g^^^^g^_^uV(l-  |-)(h/Aa)Q*] 


-  sin'*(JT{J^[i|(Aa/h)^(l-  |-)(g^f)^]ig?d?}  dr 


dg^d^i-dr 
2 


Then   substituting   <5W   into  equation   (22)   again  and 
integrating,  one  obtains 

2 

g""+  |g'"-  -^g"+  -^9'-  ^3  -   9(i-  -|-)|(g'f)' 

2    * 
=  (1-  -J-)-^  (23a) 


2  2 

in  which  A!'=(Aa/h)    and   X=/8u   are  additional  parameters 

related   to   the  amplitude   and   the  angular    frequency, 
respectively. 

Substitution  of  equations  (20a),  (20b)  and  (20c)   into 
the  governing  equation  (18b)  yields 


f..+  if._  -L-f    +  _|_(g.)2  ^  0  (23b) 
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where   the  prime  as  the  superscript  denotes  differentiation 

with   respect   to   ^,   i.e.,   f'=   df/d^.   Based   on    this 

operation,  the  time  variable  is  eliminated  from  the 
governing  differential  equations  with  an  average  minimum 
error  over  a  complete  cycle  of  the  motion.  The  assumed 
harmonic  motions  are  governed  by  the  pair  of  nonlinear 
ordinary  nonlinear  differential  equations  (23a)  and  {23b). 

In  the  same  way,  substituting  equations  (20a), (20b) 
and  (20c)  into  equation  (19),  the  boundary  conditions  can 
be  converted  to  be 


(  g  )^=i  =  0 


(  g')^.,  =  0 


(  ^'\=R  =   ° 


(24a) 


''-r  ^^=1  =  ° 


(  g"'+  |g'  '-  |yxg  )^^^   =  0 


In  order  to  obtain  a  unique  relationship  between 
amplitude,  k,  and  frequency,  X,  an  additional  restriction, 
which  supplements  the  boundary  conditions,  must  be  set. 
This  requirement  is  fulfilled  by  introducing  a 
norminalization  condition 

(  g  )^^^   =  1  (24b) 
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Finally,  The  equations  (23a),  (23b)  and  boundary 
conditions  (24a),  (24b)  which  constitute  a  nonlinear 
eigenvalue  problem  describing  the  finite  amplitude  harmonic 
response  of  an  orthotropic  circular  plate,  are  expressed  in 
terms  of  shape  functions  f,  g  and  their  derivatives. 
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Chapter  IV    NUMERICAL  ANALYSIS 

Boundary  value  problems,  as  the  name  implies,  have  the 
property  that  conditions  are  specified  at  two  sets  of 
values  of  the  independent  variables.  In  constrast  to  the 
boundary  value  problem,  the  initial  value  problems  are 
specified  by  one  set  of  values  of  the  independent  variables 
at  the  initial  point.  Consequently,  the  theory  of  the 
existance  and  uniqueness  of  solutions  to  boundary  value 
problem  is  in  a  less  satisfactory  condition  than  the 
corresponding  theory  for  initial  value  problem.  It  should 
be  expected  that  numerical  solution  of  a  boundary  value 
problem  for  a  given  ordinary  differential  equation  will  in 
general  be  a  more  difficult  matter  than  the  numerical 
solution  of  the  corresponding  initial  value  problem.  [13] 

The  nonlinear  eigenvalue  problems  formed  by  nonlinear 
ordinary  differential  equations  (23a)  and  (23b),  associated 
with  the  boundary  conditions  (24a)  and  (24b),  however,  are 
complicated.  Hence  it  is  more  convenient  to  solve  the 
problem  through  the  application  of  numerical  integration  to 
the  associated  initial  value  problem. 

For  computational  purposes,  the  associated  initial 
value  problem  is  started  by  converting  the  nth-order 
differential  equation  into  n  first  order  differential 
equations.    Then,   apply   the  Runge-Kutta-Gill  integration 
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method  [14]  to  integrate  the  differential  equation  from 
initial  point,  with  the  guessed  missing  and  the  given 
initial  values,  to  final  point.  The  value  obtained  at  the 
final  point,  in  general,  can  not  satisfy  the  given 
conditions.  This  solution  is  then  used  to  adjust  the 
initial  missing  values  by  the  Newton-Raphson  scheme  in  an 
attempt  to  make  the  next  solution  come  "closer"  to 
satisfying  all  of  the  necessary  conditions.  If  these  steps 
are  repeated  and  the  iterative  procedure  converges,  the 
boundary  conditions  (24a)  and  (24b)  will  eventually  be 
satisfied . 

a^ General  formulation . 

Any  nth-order  differential  equation,  linear  or 
non-linear,  may  be  reduced  to  a  set  of  n-simultaneous 
first-order  differential  equations  [15].  In  the  same  way, 
Equations  (23a)  and  (23b)  can  be  reformulated  to  an 
equivalent  system  of  six  first  order  equations. 


dY    -    -        * 
2|  =  H(?,Y;X,Ar,y,Q  ) 


R  <  ?  <  1 


(25a) 


where 


Y(^)  = 


■  g 

^1 

g' 

^2 

g' ' 

— 

^3 

g.  , , 

V4 

f 
f ' 

k\ 

(25b) 
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Again,    the   prime 


as 


the    superscript    donates 


differentiation,   i.e.,   g'=dg/d?-,  and  H  is  the  appropriate 
(6x1)  vector  function  defined  as  : 


H  = 


2 

-  Iy.+^Y.-  ^Yo  +  XY,+9(1-  £-)J(y  v^+y^yj 


^M   2-'3    3^2 


$    '^'^2'e    ■'3^5 


2    * 


-  -V  +  —V  -  ^(y  ) 

^^6    2^5   2^^^2' 


(26) 


The  boundary  conditions  (24a)   and   the  normalization 
condition,   y. (R)   =   1,   can  be  expressed  in  the  following 

forms  : 


■  g 

•  1 

g' 

0 

[M]Y(R)    = 

g" 

..ig.. 

-  fyXg 

"• 

0 

.  f- 

f^ 

■ 

?=R 

.    0 

(27a) 


and 


[N]Y(1)  = 


g 
g' 


?=i 


0 
0 
0 


(27b) 
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where  [M]  and  [N]  are  the  matrices 


and 


[M]  = 


1 
0 
-RXy/2 
0 


[N]  = 


1 
0 
0 


0 

1 

0 
0 


0 

1 

0 


0 
0 

1/R 
0 


0 
0 
0 


0 

0 
0 


0 

0 

1 

0 


0 

0 

-a 


0 
0 

0 

-a/R 


0 

0 

1  J 


0 
0 
0 

1 


(28a) 


?=R 


(28b) 


^  =  1 


The  first  row  of  [M]  is  the  normalized  condition, 
y  (R)=l,   at   point   ?=R  which  is  the  Junction  of  the  plate 

and  rigid  mass.  The  physical  meanings  of  the  remaining 
rows  of  [M],  referring  to  the  boundary  conditions  (24a), 
are  the  conditions  of  zero  slope,  the  shear  force  at  the 
inner  boundary  caused  by  inertial  force  of  the  rigid  mass 
and  zero  radial  displacement.  Likewise,  the  rows  of  [N] 
define  the  boundary  conditions  of  zero  transverse 
displacement,  slope,  and  radial  displacement  at  the  outer 
boundary  of  the  circular  plate,  ^=1. 

t> •  Initial-Value  Problem 

To  obtain  a  unique  solution  of  equations  (25a),  six 
boundary  conditions  must  be  specified  at  the  boundary 
points  ^=R  and  ?=1 .  However,  for  this  boundary  value 
problem  six  given   conditions   are   not   specified   at   one 
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boundary  point.  In  order  to  integrate  numerically  the  six 
first-order  nonlinear  differential  equations  by  computer, 
the  missing  conditions  must  be  specified  with  certain 
guessing  values.  These  values  are  called  missing  initial 
values.  If  the  unknown  conditions  at  one  boundary  ^=R  are 
assigned  for  certain  values,  step  by  step  integration  of 
equation   {25a)   from   ^=R   to   ^=1   would  become  possible. 

Obviously,  the  values  of  Y(^)  at  ^=1  thus  obtained  would 
not,  in  general,  satisfy  the  given  conditions  at  the  other 
end.  The  problem  is  therefore  to  determine  the  correct 
Initial  values  at  ^=R  so  that  all  the  given  conditions  at 
^=1  can  be  satisfied.  Thus,  the  corresponding  initial- 
value  problem  may  be  expressed  as 


dZ   —   -        * 

||  =  H(?,Z;X,Af,y,Q  )      R  <  ?  <  1 


(29a) 


where 


r  Z. 


Z(R)  = 


?=R 


1 
0 

^1 
RXy/2  -  )7^/R 


arjg/R 


?=R 


(29b) 


The  Z.(R),  initial-value  vector  constructed  from  the 
boundary  and  normalization  conditions  at  ?=R,  is  identical 
with  Y^{R)   in   equation   (25b),   where   i=l,2,...6,.   The 
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parameters  r]    ,ri^   and  X  in  equation  (29b)  are  unknown  miss- 


ing initial  values  and  the  eigenvalue.   Substitution  of  the 

initial   values,   Z(R),   and   equation   (28a)  into  equation 
(27a)  yields  a  system  of  four  equations. 


[M]Z(R)  = 


r  1 

0 
0 
L  0 


(30) 


?=R 


which  satisfy  the  boundary  conditions  at  ?=R. 

A   solution   of   the  initial  value  problem  thus  formu- 
lated by  equation  (29)  is  symbolically  denoted  by 

Z(^)=  Z(R)+J^  H(?,Z;n^,)72,X,*',y,Q*)d? 

=  Z(R)+J|  K(^,Z;n,K,y,Q*)d^  (31a) 

where 


f  n^   1 


17  = 


(31b) 


is   an  unknown  vector  to  be  determined.   For  any  set  of  the 

* 
given  parameters  k,    y  and  Q  ,  the  corresponding   values   rj 

,T\        and   X  can  be  solved  through  the  initial-value  problem 

such  that  the  solution  of  equation  (29)  satisfies  the  three 
given  boundary  conditions  (27b) 


[N]Z{l;n,K)    = 


r  0  1 

0 
0 


(32) 


?  =  1 
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Assume  that  Z  ( 1 ;  n  ,  at  )  is  continuously  dif  f  erentiable 
with  respect  to  17  and  k.  Then,  by  a  theorem  in  matrix 
theory,  for  the  system  of  equations  (32)  to  have  a  unique 
solution,   a  necessary  and  sufficient  condition  is  that  the 

determinant  of  the  Jacobian  matrix,  J  =  —^      •[[N]Z(l;rj,(^)|, 


is  not  zero. 


I.e., 


detl  [N]  ^  Z{1  ■,T\,K)\    *   0 

an 


Hence,  for  any  given  values  of  a-,  y  and  Q 


Y{?)=  Z(^;rj*,/^,y,Q*) 


symbolically  denotes  a  solution  to  the  eigenvalue  problems, 
where 


_* 

n  = 


*    ^ 

1 

* 

2 

* 


(33) 


Thus,  it  is  seen  that  solving  the  boundary  value 
problem  (23)  and  (24)  is  equilvalent  to  obtaining  a 
continuous  set  of  solutions  to  the  related  initial  value 
problem  (29)  which  satisfy  equation  (32). 


_  * 


A  root  r}      may   be   found   by   applying  Newton-Raphson 
iteration  method   [16]   to  solve   equation  (32).   Starting 
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with  an  initial  guess  n='?Q  and  given  parameters  k,      y,   and 

* 
Q  ,  the  iterative  covergent  sequence 


^k4.1=  ^k-^^'^k 


k=0,l,2,3, n     (34) 


is   generated,   where   Arj,  is  the  correcting  vector  and  the 
subscript  denotes  k-th  time  of  iteration. 

Neglecting   higher-order   terras  in  the  Taylor's  series 
expansion  about  rj.  ,  one  obtains 

thus 

^n^   =  -  {[N]  -|-Z(l;nij.,'f)}"^-[N]Z(l;nj^,*f) 
^^k 

=  -  {[N](J^)  J"^-[N]Z(l;njj.,A^)  (35) 

where,   at  the  k-th  step,  the  (6x3)  Jacobian  matrix  (J.)  is 

defined  as 
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(^l^k 


az 
a^ 


?=i 


az. 


az. 


az. 


Bn, 

3.2 

ax 

az^ 

3Z2 

az2 

an, 

3.2 

ax 

^^3 

^^3 

^^3 

3., 

3.2 

ax 

3Z, 

3Z, 

3Z4 

a.. 


aZ; 

a^ 


!5 
an. 


3., 


aZj 

3^ 


a.. 


9X 
3Zj 

ax 

az^ 
ax 


(36) 


?=1 


which  represents  a  change  of  final  values  with  respect  to  a 
change  in  the  initial  value  rj,  while  [N]Z{  1 ; ..  , /^) 
represents  the  k-th  error  vector  expressed  as 


[N]Z(l;nj^,A^)  = 


V^^gj^^^ 


(37) 


Therefore,  the  correcting  vector  can  be  expresseci  in  matrix 
form  : 
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tn  = 


"Arj. 

An, 
Lax 


3^1 

3^2 

ax 

az, 
3n, 

32^ 
3.2 

3Z2 

ax 

^^3 
3^1 

3^2 

3r,2 

^^3 

ax 

1 

0 

0 

0 

0 

0- 

0 
0 

1 

0 

0 
0 

0 
0 

0 
-a 

0 

1. 

3^4 

ax 

^^5 

^^5 

^^5 

Bn, 

3.2 

ax 

^^6 

^^6 

^^6 

3Z. 


3Z. 


az. 


3r?. 


3., 


ax 


-1 


J?=i. 


ri     0    0    0     0    0 

0       10      0      0      0 
0      0      0      0   -a      1 


?=1 


I.e., 


A.. 


An, 


AX 


az. 


az, 

3.^ 


3.. 


3Z, 


-a 


3Z, 


3.,       3r?, 


3Z. 
3.^ 


az, 

3^ 


-a 


azj 
3n, 


ax 


az. 

ax" 


az, 
aT 


-a 


az^ 
ax 


-1 


?=1L 


^e-^^s 


(38) 


?=1 
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which  provides  the  linear  correction  of  the  estimated  data 

ri    ,      77   and  X.   If  the  initial  estimate  of  t\    is  chosen  in  a 

_* 
sufficiently   small   neighborhood   of   the   root   17  ,    the 

_* 
sequence  (34)  will  converge  to  the  root  17  .   The  matrix  (J, 

),  must  be  evaluated  at  each  step  of  the  iteration   process 

in   order  to  generate  the  correction  vector  A77,  and  then  to 

find   the   solutions.    Thus,   the   associated   variational 
equations  are  introduced  : 


i^(i^)  =  3^(ff)  =  (^)(|^)  (39) 

d?^ax  '   3)7  ^d^^   ^-='^ax'  ax  * 

The  specific  expressions,  which  result  from  performing  the 
operations  indicated  in  equation  (39)  can  be  presented  as 
follows : 


d_,!!i,      !!2 
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3Z 

9n^ 


d?^9r? 


8Z,  ^  az. 


^'dri 


9Z  3Z  9Z 


2  az  3Z  az  az 


(40a) 


HI      3Z, 


)    = 


3r}. 


.       8Z.  ,     3Z_       - 


— >-f^"2)(37r 


d      ^^1 
df(3n7) 


8Z, 

97r 


ci?^9n2 


9Z, 

97r 


d      ^^3 
d?(9T^) 


az 

97?^ 


,       9Z,  „    9Z^       „       9Z^       ^       3Z„ 


9Z. 


d^'dn. 


e^9rj 


^2^9,2 


^3^9n2 


9rj 


•) 


2 

9Z. 


2  9Z  9Z  9Z  8Z 


(40b) 


ri         9Z^ 

d?^3no 


rl         9Z- 

— ( — ^' 


az. 


)   = 


9?? 


2 

8Z, 


=  _i(lIi)+^( ^)-^(z  w ^^ 

?^3rj2^^2*3r}2)    ^^'^2)^8^2^ 
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H   3Z,      9Z, 


d^^ax  '    9X 


ri  92„    az 
d^^ax  '    ax 


H  9Z. 
d^^ax 


az     az     az    az 

=  -f(a>r'^^(a)r)-^(a)r)^>^<a)r)^2i 

2     az     az     az     az 
^^(^-l-^fK^ax^^^^s^a^r^^^a^air^-^^e^ax^M 


d   9^5      ^^ 

df^ax  ^  ^  ax 


6 


9Z, 


az      az     az       .„ 
d^^ax  '    e^ax  '  2^ax  ^  ^^^2^^ax  ' 


(40c) 


Also,  variational  equations  at  ?=R  may  be  obtained  by 
taking  derivatives  of  Z  with  respect  to  rj,  ,  rj-  and  X 
seperately,  namely 


r9z_i 

Larj^J?=i 


0  ■ 

0 

1 

-1/R 
0 
0  , 


{41a) 


35 


0 
0 
0 
0 

1 
a/R. 


41b) 


LaxJ  ?= 


e=R 


0 
0 
0 

yR/2 

0 
0 


(41c) 


Equations  (40a),  (40b)  and  (40c),  18  first-order  equations, 
and  equations  (41a),  (41b)  and  (41c),  the   initial   values, 

constitute   the  variational  problem.   For  a  given  vector  rj 

and  given  values  of  k      and  y  ,  this  derived  problem  and  the 

initial   value  problem   (29)  may  be  numerically  integrated 

simultaneously  over  the  interval  from  ^=R   to   ?=1   with  a 

Runge-Kutta-Gill    method.    Evaluation   of   the   resulting 

solution  to  the  variational  problems  at   ^=1   provides   the 

values   for   each  entry  of   the   Jacobian  matrix   (J,), 

11 

corresponding  to  the  given  values  of  rj-  '^   and  y  .    Recall 
equation  (38),  which  may  be  rewritten  symbolically  as 

[aT,^]    =    (J^)"^[Z]  . 

Premultiply  (J^)  on  each  side  of  the  equation,  we  obtain 


(J^)  [Arjj^]  =  [Z] 
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Then,   applying  Gaussian  Elimination  method,  we  can 
solve  the  matric  equation  easily,   and   obtain  a   set   of 

correcting   vectors   Arj,   for   this  step  k.  Substitution  of 

this  correcting  vector  into  equation  (34),  we  obtain  a  new 

set   of   trial  values  for  the  vector  rj.   Repeating  the  same 

procedure  until  the  results  converges   within   a   specified 

error  bound  which  is  usually  related  to  the  accuracy  of  the 

integration  method  used. 

-* ( i  )  ( i ) 

Having   obtained  a  root  rj   ■*   corresponding  to  k=k^-', 

with  a  fixed  y,  the  solutions  corresponding  to  the  value  of 

amplitude  parameter  k^^'    are  obtained.   Next,  in  order  to 

obtain  a  solution  for  a  higher  value  of  k      than  at  "^  ,   the 
method  of   continuation  is  applied.  The  value  of  amplitude 

parameter  k^"^'    is  perturbed  to  be 

^(5  +  1)=,  ^(3)+  ^^(J)^      j=0,l,2,3,...n  . 
For  every  new  k,    iteration  is  restarted  from  the  previously 
obtained  values   of  ri=T\    ^ -^ '  .      If   Ak^-^'    is  small  enough  for 

n    "^      to  be  within  the  new  contraction  domain  of  the  Newton- 

Raphson  method,   iteration  will   converge  to  the  new  root 

-  *  (  i  + 1 )  (  i  + 1 ) 

T\    '  •"         ,    corresponding  to  the  new  k^  "^      '  .      Successive   this 

repetition  until  k   reaches  the  reasonable  large  value. 
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Chapter  V   Numerical  Computation 

The  application  of  a  numerical  integration  technique 
is  suggested  for  the  analysis  presented  previously.  Thus, 
by  integrating  the  initial-value  problem  (29)  with  a 
fourth-order  Runge-Kutta-Gill  method,  and  performing  the 
successive  iterations  by  Newton-Raphson  method, 
approximate  numerical  solutions  to  the  boundary  value 
problem,  (23)  and  (24),  can  be  obtained.  A  FORTRAN  77 
computer  program  listed  in  appendix  B  for  numerical 
calculation  is  available  for  this  particular  problem. 

The  following  procedure  of  numerical  computation  is 
suggested  and  used  in  the  investegation  of  the  problem  : 

First,   the  problem  of   free   linear  vibration   is 

* 
considered   .   Thus,  the  loading  parameter  Q  and  amplitude 

parameter  k   are  set  equal  to   zero;   the  elastic   constant 

ratios,  a  and  /9,  for  the  orthotropic  plate  are  set  equal  to 

a  specific  value.   In  this  case  the  equation   (23a),   which 

governs   the  transverse  displacement,  becomes  linear.   Then 

the   initial-value   problem    (29)    and    the    associated 

variational   problem   (40)   and   (41)    can   be   integrated 

numerically  over  the  interval  [R,l]  with  a  set  of  estimated 

initial    values    r?^ .     By    the   Newton-Raphson  method, 
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successive  correction  and  integration  are  carried  out  until 
the  error  norm  satisfies  the  inequality  : 

Max| [N]Z{1) I  <  O.lxio"^ 

where  Maxi  i ,  the  maximum-element  norm  of  the  error  vector, 
is  consistent  with  the  forth-order  Runge-Kutta-Gill 
integration  method  with  step  size  /j=l/40.  Therefore,  the 
solution  for  the  linear   free   vibration   of   the   circular 

_* 
plate   is   obtained  and  the  corresponding  values  of  17   are 

stored.   From  this  solution  of  the   linear   free   vibration 

the  solutions  of  finite  amplitudes  can  be  examined  with  the 

concept  of   neighboring  solutions,   which   is   called   the 

method   of   continuation.    The  discrete  representations  of 

the  resonance  curves  and  accompanying  solutions   are   found 

by  successively   increasing  the  value  of  amplitude,  *■ ,  and 

re-starting  the  correction  and  integration  procedure   from 

those  values   of   (r?^,  n    ,      X)   obtained   in  the  solution 

corresponding  to  the  previous  value  of  k.      This  process   is 

terminated  when  k   reaches  a  reasonable  large  value,  k^ ,    and 

the  corresponding   initial   values   of   (17^,  ru>      ^"')   are 

stored.  At  this  stage  the  steady  state  response  due  to 
forced  vibration   can  be   determined  by  a   perturbation 

technique  .  On  setting  k=k^ ,  {r)^,  n^.  X)  =  (rj^,  xj^,  x"*)  , 
and  applying   a   small   load,   Q  (^),   to   the   plate,   the 
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response  is  easily  determined  with  a  slight  modification  of 
the  above  mentioned  process.   Namely,  the  type  of   loading, 

« 
Q  {^)   is   held   fixed  while   the  value  of  k    is  gradually 

decremented  from  k  as  the  integration  and  correction 
procedure  of  the  initial-value  method  yields  successive 
solutions  to  the  forced  problem.   Hence,  distinct   response 

curves  are  obtained  by  merely  changing  the  sign  of  Q  {?). 
The  totality  of  this  approach  is  illustrated  schematically 
in  figure  2. 
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Figures  3  through  6  represent  the  graphs  of  the 
dimensionless  angular  frequency  vs  the  elastic  constant 
ratio,  c.  The  results  point  out  that  the  frequency 
parameter  cj  for  an  axisymmetric  free  vibration  of  an 
orthotropic  circular  plate  carrying  a  concentric  mass 
increases  as  the  stiffness  parameter,  c,  increases,  while 
the  radius  ratio  is  held  constant.  These  data  provide  the 
starting  frequency  of  nonlinear  free  vibration  for 
different  materials  which  possess  different  stiffness 
parameters.  If  the  stiffness  constant  becomes  unity,  that 
is  c=1.0,  then  the  orthotropic  plate  case  is  simplied  to  a 
isotropic  plate  case.  The  results  of  linear  free  vibration 
of  isotropic  circular  plate  with  different  mass  ratio  are 
illustrated  in  Table  1.  Close  agreement  with  the  results 
found  by  Becker  [5]  and  Chiang  and  Chen  [4]  is  observed 
directly. 

The  results  in  tables  2  through  5  list  some  of  the 
data  set  in  figures  3  through  6,  which  illustrate  the 
influence  of  variations  of  stiffness  parameters  c  on  linear 
angular  frequency.  The  higher  stiffness  parameter  causes  a 
higher  angular  frequency  as  the  mass  ratio  is  held  fixed. 

Figures  7  and  8  present  the  results  determined  for 
nonlinear  free  vibration  with  various  values  of  stiffness 
parameter  ( c=0 . 5 ,  1.0,  and  2.0),  radius  ratio  ( R=0 . 1  and 
0.3).   At  the  same  amplitude,  the  materials  with   a   higher 
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stiffness  parameter  oscillate  at  higher  resonant  frequency, 
while  the  materials  with  a  low  stiffness  parameter 
oscillate  at  a  lower  resonant  frequency.  In  addition,  the 
lower  stiffness  material  has  a  steeper  resonant  curve  than 
higher  stiffness  materials.  The  physical  meaning  of  this 
phenomenon  is  that  the  resonant  frequency  of  higher 
stiffness  materials  increases  faster  than  that  of  lower 
stiffness  materials  while  the  amplitude  increases  at  the 
same  amount . 

In  figures  9-14,  forced  nonlinear  vibrations 
areindicated  with  various  values  of  stiffness  parameters 
(c=0.5,  1.0,  and  2.0),  radius  ratios  {R=0.1   and   0.3)   and 

* 
acting   force   parameters   (Q  =0.0,   -20.0,  20.0,  -40.0  and 

40.0).   In  practice  it  is  the  custom  to  plot  iw/hi   against 

u   .    It   is  noted  also  that   the  quantity  of  amplitude  is 

negative  on  the  response  curves  to  the  right   of   the   free 

* 
oscillation  curve  for  Q  =0  and  positive  to  the  left  of  it, 

which  means  that  the  motion  is  in  phase  with   the   external 

*  • 

force  when  Q   is  postiive  and  180   out  of  phase  with  the 

external  force  when  Q   is  negative.  This   behavior   of 

nonlinear   resonance   is  similar  to  that  found  for  the  free 

and  forced  oscillation  of  Duffing 's  hard  spring  system   for 
a  single  degree  of  freedom  [17]. 
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Chapter  VI    Stress  Analysis 

Stress  distribution  on  the  circular  plate  can  be 
obtained  from  the  stress-strain  relation  and  strain- 
displacement  relation  [7].  The  following  symbols  and 
expressions  are  used  to  denote  the  stresses 


o    :  radial  bending  stress 
r  ^ 


a  -:    circumferential  bending  stress 


a"*:  radial  membrane  stress 

r 


e 
where 


o^:    circumferential  membrane  stress 


b  /Sh      ,     ^  /8    , 

®      2a22(/8-a2)a    ^^   ^   ^ 


a'"=  +  -^— (I) 


r   -   a22  ^? 


and       o   =   +  — (4>^) 


which  are  derived  from  equations  (1)  and  (2),  relate  the 
bending  and  membrane  stresses  to  the  dimensionless 
deflection,  X'  and  stress  function,  <i> ,    respectively. 
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In  terms  of  the  previous  assumptions. 


X(?,T)=  A  g(^)sinoT 


2        2 

<t>{^,T}=   A  f(^)sin  or 


and  when  time,  t,  is  equal  to  an  odd  multiple  of  n/2u      the 
maximum  excursions  occur,  thus  stresses  become 


b    2 

<y   a„„a 


r^22^       /8  /"aT   ,  ,  ,  ,  a  ,  , 
h       "   2(/S-a^)        ^ 


^b^   ^2 

Ie!|2L  =  +  _^_Zl__(crg..+  fgM  (44) 

h'^       ~  2(fi-a)  ^ 


m    2 
^^^22^    ^   ,f, 


2       -   ^? 


m    2 
and       5 ~   t   '^^^ 


Figures  15  through  17  represent  the  radial  bending 
stresses;  figures  18  through  20  represent  the  radial 
membrane  stresses;  figures  21  through  23  represent  the 
circumferential  bending  stresses,  while  figures  24  through 
26  represent  the  circumferential  membrane  stresses.  All  of 
these   stresses   are   calculated  with  radius  ratio,  R=0 . 3 , 
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mass  ratio,  y=1.0,  in   the   condition  of   free  vibration, 

i.e.,  Q  =0.0.  Also,  the  stress  distribution  at  the  inner 
edge  and  at  the  outer  edge  of  radial  membrane  stress  and 
radial  bending  stress  are  shovm  in  figures  27  through  30 

From  equations  (44),  it  is  noted  the  amplitude  has  a 
great  influence  upon  the  distrubution  of  bending  stress. 
The  system  behaves  as  a  linear  property  when  the  amplitude 
is  small,  say,  /c=0.2.  On  the  contrary,  a  higher  amplitude 
causes  stress  to  be  larger  and  has  a  nonlinear 
distribution.  Also,  elastic  constant  ratios  a  and  ^,  shape 
functions  f  and  g  and  their  derivatives  and  location  ?, 
play  an  important  role  on  the  stress  distribution. 
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Table  2 
Dimensionless  Angular  Frequency  u 


y  =  1.0 

0.1 

0.2 

0.3 

0.4 

0.5 

8.1811 

8.5052 

9.2413 

10.5576 

1  .0 

10.4410 

11.1341 

12.3858 

14.4140 

1.5 

12.2449 

13.2373 

14.8760 

17.4362 

2.0 

13.8016 

15.0459 

17.0044 

20.0067 

2.5 

15.1942 

16.6580 

18.8941 

22.2826 

3.0 

16.4667 

18.1267 

20.6112 

24.3465 

Table  3 
Dimensionless  Angular  Frequency  o 


r  =  2 .0 

0.1 

0.2 

0.3 

0.4 

0.5 

8.0053 

7.9246 

8.1674 

8.8931 

1.0 

10.1994 

10.3559 

10.9334 

12.1348 

1.5 

11.9530 

12.3031 

13.1261 

14.6764 

2.0 

13.4673 

13.9787 

15.0009 

16.8385 

2.5 

14.8224 

15.4728 

16.6659 

18.7528 

3.0 

16.0610 

16.8343 

18.1789 

20.4891 
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Table  4 
Dimensionless  Angular  Frequency  cj 


y  =  8.0 

0.1 

0.2 

0.3 

0.4 

0.5 

7.1239 

5.8871 

5.3415 

5.3524 

1  .0 

9.0084 

7.6526 

7.1365 

7.2987 

1.5 

10.5241 

9.0746 

8.5618 

8.8254 

2.0 

11.8730 

10.3007 

9.7813 

10.1243 

2.5 

13.0139 

11 .3950 

10.8646 

11.2746 

3.0 

14.0908 

12.3928 

11.8493 

12.3179 

Table  5 
Dimensionless  Angular  Frequency  u 


y  =  16.0 

0.1 

0.2 

0.3 

0.4 

0.5 

6.2795 

4.6424 

4.0171 

3.9324 

1.0 

7.8941 

6.0229 

5.3642 

5.3615 

1.5 

9.2006 

7.1371 

6.4345 

6.4827 

2.0 

10.3351 

8.0985 

7.3503 

7.4367 

2.5 

11.3536 

8.9569 

8.1639 

8.2815 

3.0 

12.2864 

9.7398 

8.9035 

9.0477 
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'  '  I  '  '  1  '  '  I  '  '  I  '  ■  I  <  '  I  I  I  1  I  I  I  I  I  I  I  I  l' 
0.0  0.3  0.6  0.9  1.2  1.5  1.8  2-1  2.4  2.7  3.0 


Fig.  3.   Dimenslonless  Angular  Frequency  with  R=0.1 
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0.0  0.3  0.6  0.9  1.2  1.5  1.8  2-1  2-4  2.7  3.0 


Fig.  4.   Dimensionless  Angular  Frequency  with  R=0 . 2 
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u. 


'  '  I  '  ■  I  ■  '  I  '  '  I  '  '  I  ■  I  I  I  ■  I  I  I  I  I  I  I  I 
0.0  0.3  0.6  0.9  1.2  1.5  1.8  2.1  2.4  2.7  3.0 


Fig.  5.   Dimensionless  Angular  Frequency  with  R=0 . 3 
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'I  I  '  I  '  '  I  ■  '  I  '  ■  1  '  '  I  ■  '  I  ■  '  I  '  '  I  ■  '  I  '  '  I 
0.0  0.3  0.6  0.9  1.2  1.5  1-8  2-1  2.4  2.7  3-0 


Fig.  6.   Dimensionless  Angular  Frequency  with  R=0 . 4 
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Fig.  15.   Radial  Bending  Stress  with  c=0 . 5 
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63 


o 

•  CO  o  n 

0  *  *  • 

1  O  t-i  o 
•  I  H  R 


M^ 


o 

N 


o 


00 

o 


o 


«0 

o 


o 


o 
II 
U 


to 

c 
« 

h 

JQ 
E 
0) 


10 


I 


0» 


I  ■  ■  '  '  1  ■  ■  •  ■  I 

0»  00  N. 


(O 


»o 


'    I    I    I    I    I    I   I    I    I    I 


o 


CM 

« 
CM 
CM 

e 
B  u 


CM 


64 


o 

•  to  O  o 

o  •  •     • 

U  O  «i  o 

•  II  II   n 


n 


CM 

o 

H 


n 
i|iiif|iiii|iiiii'  -^ 
r)  CM  —  o 


JLU 


o 


00 

o 


o 


o 


in 

o 


It 
o 

Xi 
4* 


CD 
« 
« 
U 

m 

« 
« 

e 
« 

s 


10 


« 
W 
CM 

(0 


CM 


65 


o 

•  pj  o  o 

o  •  •  • 

II  O  rH  O 

t  II  II  II 


JU/> 


CM 

O 

II 


at 
o 


CO 

o 


o 


(£> 

o 


o 


CM 

n 
u 


« 

(4 

0) 
01 

c 

e 
1} 

s: 


•o 
as 


o 

CM 


•H 
fa 


o 


CM 

nJ 

CM 
C>4 

S    ^ 


66 


b     2 


13  - 

• 

» 

1  1  - 

Q    =0.0 
iv=0.3 

1 

' 

y=1.0 

/ 

9- 

R=0.3 

/ 

7  - 

/ 

5- 

/ 

3- 

/ 

1  - 

J. 

>^ 

-1  - 

^,^^<^/ 

^---^        /f 

=0 

.2 

■^^""'^     y^     / 

-3- 

— 

•                            / 

-5- 

=1 

.0 

,^ 

■          / 

-7  - 

/ 

1 
-9- 

/ 

1  1  - 

Ar= 

4. 

0  ^y^ 

13  -^ 

1      1      1      1      1 

— I — 

T — r 

1 1 1 1 1 1 1 

............. 

1   1   I   1   I   1 

0.3    0.4    0.5    0.6    0.7"    0-8    0-9    1.0 
Fig.  21.   Circumferential  Bending  Stress  with  c=0 . 5 


67 


13 


1  1 


b    2 


7  - 
5- 
3- 


1  1  - 


-13- 


Q  =0.0 
1^=0.3 
y=1.0 
R=0.3 


I  I  I  I — I — I — I — I — [— I — I — r— I — I — I— I — I — I — I — I — I — ]  I  [ — 1—1 — I — I  I  I  I  I  I 

0.3    0.4    0.5    0.6    0.7    0.8    0.9    1.0 


Fig.  22.   Circumferential  Bending  Stress  with  c=1.0 


68 


13 


1  1 


b     2 
<j  a   a 
e  22 


9 

5 
3 

1 


-1  - 


-3- 


-5- 


-7  • 


-1  1  - 


13- 


/f=1.0 


K=A.O 


Q  =0.0 
w=0.3 
y=1.0 
R=0.3 


I  I  I — I — [  I  -I — 11  I 


I  I  r 


~i — T""' — ' — ' — ' — I — ' — ' — ' — ' — I — ' — ' — I""" — I — I — '""' — I — r 


0.3    0.4    0.5    0.6    0.7    0.8    0.9    1.0 
Fig.  23.   Circumferential  Bending  Stress  with  c=2 . 0 


69 


«/• 


o 

•  m  o  n 

o    •  •    • 

II  o  •-•  o 

«      II  II   II 

Of  :i  >•« 


a> 
o 


00 

o 


o 


o 


o 


o 


o 


lO 


(0 
(0 
ID 

0) 
« 

c 

(4 

6 
0) 


V 

u 


e 

u 
u 


til 

•ri 
{ki 


CM 

(0 

CM 

CM 

CM 

(0 

JC 

e  <x> 

to 

70 


o 

'  n  o  ci 
o    •    •    • 

II      O    r-l    O 
»  II       II       II 

a  i  >•« 


II 


<o 


r— 1 1 r- 


JUl^ 


cr> 

iH 

• 

II 

o 

U 

Xi 

*J 

■H 

% 

00 

0) 

o 

(0 

O 


U3 

o 


IT) 
O 


^4 
■♦J 
10 

c 

(0 

£i 

B 
O 

s 


(0 

•H 

h 
V 

H 

O 
U 


Mi 

CM 


0) 

■H 


o 


CM 

(D 

CM 
CM 

(0 

6  a> 
b 


71 


o 

•  PJ  O  CO 

o  •  •  • 

U  O  f-t  o 

»  II  II  II 


H/- 


o 


00 

o 


O 


o 


in 
o 


II 


CM 

o 
II 


M 


Xi 

n 
o> 
tt 

CO 

c 

(0 

XI 

6 


ta 

c 

91 
U 
V 
<M 

6 

o 

h 
O 


(0 

CM 


Ok 


I    '     '    '     ■     I    '     '    '     ■    I     I     I    I     I    I     I     I     I     I    I 

1^  to  m  •v  r> 


<v 


O 


n 
o 


(0 
CM 
CM 

(0 

O 


CM 

J2 


72 


^>22a^ 

h^ 

10- 

Q    =0.0 

, 

^/=0.3 

9- 

y=1.0 

. 

R=0.3 

-  c=2.0 


c=1.0 


1  '  I  ' — I — ' — I — ' — I — ' — I — ' — I — > — I — ' — \ — I — I — I — r 
0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0 


Fig.  27,   Radial  Membrane  Stress  at  the  Inner  edge 
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CONCLUSION 

The  energy  method  for  deriving  basic  governing 
differential  equations  has  been  presented  in  this  thesis. 
The  merit  of  this  approach  is  that  it  provides  automatically 
the  governing  differential  equations  and  the  correct  number 
of  boundary  conditions  and  their  expressions. 

An  exact  solution  for  the  non-linear  boundary-value 
problem  is  very  complicated.  Thus,  a  approximate  method 
must  be  used.  Based  on  the  assumption  of  the  existence  of 
harmonic  oscillations,  the  Ritz-Kantorovich  averaging  method 
is  introduced  to  convert  the  partial  differential  equations 
to  ordinary  differential  equations  by  elimination  of  the 
time  variable,  which  makes  the  problem  easy  to  solve.  Then, 
approximate  numerical  results  can  be  obtained  by  the 
initial-value  problem  method,  the  related  Runge-Kutta-Gill 
integration  method,  and  the  Newton-Raphson  iteration  method. 
The  use  of  numerical  integration  leads  to  a  discrete 
approximation  of  the  continuous  system. 

It  has  been  noted  that  the  elastic  constant  ratio,  c, 
and  mass  ratio,  ,  have  great  influence  on  the  vibration. 
At  the  same  value  of  amplitude,  the  higher  frequency 
accompanies  the  higher  stiffness  and  lower  mass  ratio  ;  on 
the  contrary,  the  lower  frequency  is  related  to  the  lower 
stiffness  and  higher  mass  ratio. 
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APPENDIX   A 

The  generalized  Hooke ' s  law  for  the  case  of 
cylinderical  orthotropy  can  be  writen  in  the  form: 

e    11  e    12  r 

where  a. .  are  elastic  constants  (compliance  coef f iciences) . 

In  general,  the  strain  in  a  direction  normal  to  an 
applied  tensile  stress  will  be  negative  and  the  change  in 
the  voliime  of  an  element  subjected  to  a  tensile  stress  will 
be  positive.  Thus,  under  tensile  stress,  the  change  of 
volume  V   in  an  element  will  obey  the  inequality 

Vo(e  +  €+    €  )  >  0  (A2) 

0^  r   e   z 

when  there  is  only  a  radial  tensile  stress  o        the  volume 

r 

change  can  be  expressed  as : 

Vr(^12-'  ^22-^  ^23^  ^  °  ^^^^ 

and  for  a  circumferential  tensile  stress  c 

Ve^^ii-'  ^12-*-  ^13^  >  °  (^^) 

equations  (A3)  and  (A4)  can  be  rewritten  as  follow: 

"''-<1  (A5) 


^22      ^22 
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^^2      ^13  ,  ,  ,  (A6: 


^11         ^11 


But,  a    a   ,  and  a    are   negative.    Thus   "^12'^^ 


22 


and 


-a,„/a,,  are  restricted  to  values  between  0  and  1;  i.e., 
12   11 

^12 
0  < ^^^  <  1  (A7) 

^22 

^12 
0  < ^^^   <  1  (A8) 

^11 
from  previous  definitions  on  page  7,  a  =  "^lo^^ll  ^^^   ^   ~ 
a   /a   ,  the  restrictions  in  equations  (A?)  and  (A8)  can  be 

rewritten  as: 

0  <  a  <  1  (A9) 

Also,  Poisson's  ratio  u   =   -   ^yo^^ll'    ^^^^^^^      °^  page 
7,  must  be  less  than  one.  Hence 

V   =  (-312/^11^/^^22/^11^  =  a//9  <  1 
or  a.    <   ^  (AlO) 
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APPENDIX  B 


C * 

c  * 

C  FREE  VIBRATION  OF  A  CLAMPED  ORTHOTROPIC  CIRCULAR  PLATE  * 

C  * 

C  WITH  CONCENTRIC  RIGID  MASS  AT  THE  CENTER  * 

C  * 

C  R=0.3   VV=0.3   GAMMA=1.0  BETA=  'VARIABLE'  * 

C  * 

C * 

C 

C  VARIABLES  DESCRIPTION: 

C 

C  A:  AMPLITUDE  PARAMETER 

C  DA:  INCREMENT  IN  AMPLITUDE 

C  R:  RADIUS  RATIO 

C  DR:  INCREMENT  IN  RADIUS  RATIO 

C  H:  STEP-SIZE  FOR  NUMERICAL  INTEGRATION 

C  W:  POISSON'S  RATIO 

C  GAMMA:  MASS  RATIO 

C  DG:  INCREMENT  IN  MASS  RATIO 

C  ALPHA:  ELASTIC  CONSTANT  RATIO 

C  BETA:  ELASTIC  CONSTANT  RATIO 

C  C:  ELASTIC  CONSTANT  RATIO  (INVERSE  OF  BETA) 

C  QSTART:  LOADING  DENSITY 

C   ***************************************************** 

IMPLICIT  DOUBLE  PRECISION   (A-H,0-Z) ,INTEGER(I-N) 
INTRINSIC  ABS.SQRT 

DIMENSION  Y(24) ,Q(24) ,TP(3,4) ,C(3) ,ER(3) ,D(6,41) 
DIMENSION  SP(50) .SRA(50) ,FRE(50) ,AMP(50) 
CHARACTER  *5  RR 

112  F0RMAT(5X,I2,3X,'FREQ='.F8.4,3X,'AMP=',F8.4) 

113  FORMAT(/,10X,'QSTAR=',F6.2,/) 

114  FORMAT(/,10X,'BETA=',F5.2,/) 

119  FORMAT(/,10X,'RADIUS  RATIO=',F6.2 ,15,/) 

120  FORMAT(/,10X,'MASS  RATIO=',F6,2,/) 
QSTAR=O.OD-0 

A«O.OD-0 

P=10.2158D-0**2 

BETA=2.0D-0 

GAMMA=O.OD-0 

ALPHA=0.3D-0*BETA 

DR=O.OD-0 

DG=1.0D-0 

W=0.3D-0 

R=0.3D-0 

H=2.5D-2 
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LL=29 

WRITE  (6,113)QSTAR 

WRITE(6,114)BETA 

DO  510   IR=1,1 

LL=LL-0 

R=R+DR 

WRITE(6,119)R,LL 

DO  520    IG=1,1 

IK=1 

GAMMA=GAMMA+DG 

WRITE(6,117) 

WRITE (6, 120) GAMMA 

WRITE(6,117) 

WRITE(6,117) 
C  ***  CONSTRUCT  INITIAL  VALUES 
500   DO  10  1=1,24 
10    Y(I)=0.0D-0 

Y(l)=1.0D-0 

Y(3)=-4.67D-0 

Y ( 4 ) =-Y ( 3 ) /  R+0 . 5D-0*R*GAMMA*P 

Y(5)=0.82D-0 

Y(6)=(ALPHA*Y(5))/R 

Y(9)=1.0D-0 

Y(10)=-1.0D-0/R 

Y(17)-1.0D-0 

Y(18)»ALPHA/R 

Y(22)=»0.5D-0*R*GAMMA  — 

IF(IK.EQ.l)  GO  TO  600 

DO  15   1=1,6 
15         Y(I)-D(I,1) 
C  ***  X=INDEPENDENT  VARIABLE 
600       X=R 

DO  20   1=1,24 

20  Q(I)=O.OD-0 
DO  21    1=1,6 

21  D(I,1)=Y(I) 

C  ***  PERFORMANCE  "RUNGE-KUTTA-GILL"  INTEGRATION 

DO  25  1=2, LL 

CALL  RKGPL(X,H,Y,Q,P,A,ALPHA,BETA,QSTAR) 

DO  30  J=l,6 
30    D(J,I)=Y(J) 
25    CONTINUE 
C  ***  ER(I)=ERROR  VECTOR  FOR  BOUNDARY  CONDITION  AT  X=1.0 

ER(1)=D(1,LL) 

ER(2)=D(2,LL) 

ER(3)=D(6,LL)-ALPHA*D(5,LL) 

DO  35  1=1,3 

DER=DABS(ER(I)) 

IF(DER.GT.0.1D-5)  GO  TO  36 
35    CONTINUE 
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GO  TO  900 
36    CONTINUE 

C  ***  TPd.J)  IS  THE  JACOB  IAN  OF  THE  MAPPING  OF  INTERAL  VALUES  TO 
C     FINAL  VALUES 

TP(1,1)=Y(7) 

TP(2,1)=Y(8) 

TP(3,1)=Y(12)-ALPHA*Y(11) 

TP(1,2)=Y(13) 

TP(2,2)=Y(14) 

TP(3,2)=Y(18)-ALPHA*Y(17) 

TP(1,3)=Y(19) 

TP(2,3)=Y(20) 

TP(3,3)=Y(24)-ALPHA*Y(23) 

DO  40  1=1,3 
40    TP(I,4)=ER(I) 

CALL  GAUSSX(TP,C,3,4) 
C  ***   C(I)=CORRECTION  VECTOR 

DO  76    1=1,6 
76  Y(I)=D(I,1) 

Y(3)=Y(3)-C(1) 

Y(5)-Y(5)-C(2) 

P=P-C(3) 

DO  80   1=7,24 
80  Y(i)=O.Od-0 

Y(4)=-Y(3)/R+0,5D-0*R*GAMMA*P 

Y(6)=(ALPHA*Y(5))/R 

Y(9)=1.0D-0 

Y(10)=-1.0D-0/R 

Y(17)=1.0D-0 

Y(18)=ALPHA/R 

Y( 22 ) =0 . 5D-0*R*GAMMA 

GO  TO  600 
900        SRA(IK)=DSQRT(A) 

PP=ABS(P)/BETA 

SP(IK)=DSQRT(PP) 

WRITE(6,112)IK,SP(IK),SRA(IK) 

IF    (A.LE.0.09D-0)THEN 

DA=0.01D-0 

ELSE 

DA=0.1D-0 

END  IF 

A=A+DA 

IK=IK+1 

IF(IK.GT.50)GO  TO  520 

GO  TO  500 
520        CONTINUE 

OPEN(UNIT=10,FILE='RR') 

REWIND   10 

DO  44    IK=1,50 

FRE(IK)=SP(IK) 
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AMP(IK)=SRA(IK) 

WRITE(10,45)FRE(IK),A>IP(IK) 
45    FORMAT(5X,2F10.4) 
44    CONTINUE 

CLOSE  (10) 
510   continue 
550   CONTINUE 

STOP 

END 


SUBROUTINE  RKGPL(X,H,Y,Q,P,AP,ALPHA,BETA,QSTAR) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) ,INTEGER(I-N) 

DIMENSION  Y(24),Q(24),DY(24),A(2) 

A(l)=. 2928932188134524 

A(2)=l. 707106781186 547 

H2=0.5D-0*H 

CALL  DERIVL(X ,H,Y ,DY ,P ,AP ,ALPHA,BETA,QSTAR) 

DO  13  1=1,24 

B=H2*DY(I)-Q(I) 

Y(I)=Y(I)+B 

Q(  I ) =Q( I ) +3 . OD-0*B -H2*DY ( 1 ) 
13  CONTINUE 

X=X+H2 

DO  60   J=l,2 

CALL  DERIVL(X,H,Y,DY,P,AP,ALPHA,BETA,QSTAR) 

DO  20   I»l,24 

B=A(J)*(H*DY(I)-Q(I)) 

Y(I)»Y(I)+B 

Q( I) =Q( I) +3 .OD-0*B-A( J)*H*DY( I) 
20    CONTINUE 
60    CONTINUE 

X=X+H2 

CALL  DERIVL(X,H,Y,DY,P,AP,ALPHA,BETA,QSTAR) 

DO  26    1=1,24 

B=.166666666666666*(H*DY(I)-2,0D-0*Q(I)) 

Y(I)=Y(I)+B 

Q(I)=Q(I)+3,0D-0*B-H2*DY(I) 
26  CONTINUE 

RETURN 

END 
c 

SUBROUTINE  DERIVL(X ,H,Y ,DY ,P ,AP , ALPHA, B ETA, QSTAR) 

IMPLICIT  DOUBLE  PRECISION   (A-H,0-Z) ,INTEGER(I-N) 

DIMENSION  Y(24),DY(24) 

Cl=l .0D-0-ALPHA**2/BETA 

DO   10    1=1,3 
10  DY(I)=Y(I+1) 

DY(5)=Y(6) 

DO  15    1=7,9 
15  DY(I)=Y(I+1) 
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DY(11)=Y(12) 

DO  20   1=13,15 
20  DY(I)=Y(I+1) 

DY(17)=Y(18) 

DO  25   1-19,21 
25  DY(I)=Y(I+1) 

DY(23)=Y(24) 

IF    (AP,EQ.O.OD-0)THEN 

DY(4)=-2.0D-0*(Y(4)/X)+BETA*Y(3)/(X*X)-Y(2)*BETA/(X*X*X)+P*Y(1) 
&+9.0D-0*AP*Cl*(Y(3)*Y(5)+Y(2)*Y(6))/X 

GO  TO  90 

ENDIF 

DY(A)=-2.0D-0*(Y(4)/X)+BETA*Y(3)/(X*X)-Y(2)*BETA/(X*X*X)+P*Y(1) 
&+9.0D-0*AP*Cl*(Y(3)*Y(5)+Y(2)*Y(6))/X+QSTAR*Cl/DSQRT(AP) 
90  DY(6)=-Y(6)/X+BETA*Y(5)/(X*X)-(Y(2)*Y(2)*BETA)/(2.0D-0*X) 

DY(10)=-2.0d-0*(Y(10)/X)+BETA*Y(9)/(X*X)-Y(8)*BETA/(X*X*X)+P*Y(7) 
&+9.0D-0*AP*Cl*(Y(5)*Y(9)+Y(3)*Y(ll)+Y(2)*Y(12)+Y(6)*Y(8))/X 

DY(12)=-Y(12)/X+BETA*Y(11)/(X*X)-(Y(2)*Y(8)*BETA)/X 

DY(16)=-2.0D-0*(Y(16)/X)+BETA*Y(15)/{X*X)-Y(14)*BETA/(X*X*X) 
&+9.0D-0*AP*Cl*(Y(3)*Y(17)+Y(5)*Y(15)+Y(2)*Y(18)+Y(6)*Y(14))/X 

&+P*Y(13) 
DY(18)=-Y(18)/X+BETA*Y(17)/(X*X)-(Y(2)*Y(14)*BETA)/X 
DY(22)=-2.0D-0*(Y(22)/X)+BETA*Y(21)/(X*X)-Y(20)*BETA/(X*X*X) 

&+Y(l)+9.0D-0*AP*Cl*(Y(3)*Y(23)-»yY5)*Y(21)+Y(2)*Y(24)4yY6)*Y(20))/X 

&+P*Y(19) 
DY(24)=-Y(24)/X+BETA*Y(23)/(X*X)-(Y(2)*Y(20)*BETA)/X 

70  RETURN 

END 
c — —- — — — «____—__-._—_————  -     -       — 

SUBROUTINE  GAUSSX(A,X,N,N1) 

IMPLICIT  DOUBLE  PRECISION   (A-H,0-Z) ,INTEGER(I-N) 

DIMENSION  A(N,N1),X(N) 

DO  200   J=1,N 

J1=J+1 

IF   (Jl.GT.N)    GO  TO  980 

BIG=DABS(A(J,J)) 

M=J 

DO   900   L=J1,N 

IF   (DABS(A(L,J)).LE.BIG)    GO  TO  900 

M=L 

BIG=DABS(A(L,J)) 
900        CONTINUE 

DO   990   JJ=J,N1 

DUMMY=A(M,JJ) 

A(M,JJ)=A(J,JJ) 
990    A(J,JJ)-DUMMY 
980   CONTINUE 

S=1.0D-0/A(J,J) 

DO  201  K=J,N1 
201    A(J,K)=A(J,K)*S 


86 


DO  202  1=1, N 

IF(l.EQ.J)  GO  TO  202 

AIJ=-A(I,J) 

DO  204  K=J,N1 
204   A(I,K)=A(I,K)+AIJ*A(J,K) 
202   CONTINUE 
200   CONTINUE 

DO  300  1=1, N 
300   X(I)=A(I,N1) 

RETURN 

END 
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ABSTRACT 


The  problem  of  free  and  forced  nonlinear  large 
amplitude  vibrations  of  a  clamped  orthotropic  circular  plate 
with  a  concentric  rigid  mass  is  studied.  Hamilton's 
principle  is  applied  to  constitute  the  governing 
differential  equations  and  boundary  conditions.  Assuming 
the  existence  of  harmonic  vibration,  the  time  variable  can 
be  eliminated  by  employing  a  Kantorovich  averaging  method. 
Then,  the  governing  differential  equations  for  the  problem 
are  reduced  from  nonlinear  partial  differential  equations  to 
non-linear  ordinary  differential  equations.  The  related 
approximate  numerical  results  can  be  obtained  by  using  the 
initial-value  problem  in  conjunction  with  the  Newton-Raphson 
iteration  scheme.  The  results  reveal  the  effects  of  finite 
amplitude  and  orthotropy  of  materials  upon  the  dynamic 
responses . 

The  characteristics  of  nonlinear  vibration  as  well  as 
radial  and  circumferential  stresses  are  studied  for  various 
mass  ratios,  radius  ratios  and  elastic  constant  ratios.  For 
the  purpose  of  comparison,  the  present  orthotropis  circular 
plate  problem  can  be  reduced  to  a  flat  isotropic  circular 
plate  problem  by  setting  both  elastic  constant  ratio  and 
mass  ratio  equal  to  1.  The  results  are  in  agreement  with 
prior  work.  [3,4] 


